\section{仿真方法}\label{method}
%The finite element method (FEM) is employed in the solution of actuator's magnetostatic problem. To make it easier to simulate and analyze, in Fig. \ref{simple2}, some simplifications are applied to the model. Given the axisymmetric nature of actuator's structure, we use one half of it for calculation, as shown in Fig. \ref{bd}.
对于接触器的静磁场求解，我们采用有限元方法进行计算。为了方便仿真计算和分析，如图 \ref{simple2} 所示，我们队模型进行了相应的简化。考虑到接触器的轴对称特性，我们对二分之一的模型进行建模计算，如图 \ref{bd}。
%\begin{figure}[ht]
%    \centering{
%
%            }
%    \caption{Model for actuator's magnetostatic problem
%    }
%    \label{magnetfield} %% label for entire figure
%\end{figure}

\subsection{有限元方法}
\subsubsection{控制方程}

%In the actuator, the magnetostatic field obeys the following equation:
在接触器中，静磁场遵循如下方程：
\begin{equation}\label{possion}
  \nabla\times\left(\frac{1}{\mu_{0}}\nabla\times\boldsymbol{A}-\boldsymbol{M}\right)=\boldsymbol{J}
\end{equation}
%where $\mu_0$ is permeability of air, $\boldsymbol{A}$ is magnetic field potential vector, $\boldsymbol{M}$ is magnetization vector, and $\boldsymbol{J}$ is current density.
其中，$\mu_0$为空气的磁导率，$\boldsymbol{A}$为磁场磁势向量，$\boldsymbol{M}$为磁化向量，$\boldsymbol{J}$为电流密度。
在柱坐标 $(r,z,\phi)$中，对于轴对称结构，我们可以得到
\begin{equation}\label{axis}
  \left\{
  \begin{aligned}
  A_r&=A_z=0,&A_\phi=&A,\\
  J_r&=J_z=0,&J_\phi=&J.
  \end{aligned}
  \right.
\end{equation}
%Then \eqref{possion} is simplified as
公式\eqref{possion}简化为
\begin{equation}\label{huanjian1}
\begin{aligned}
  &\frac{\partial}{\partial z} \left[ \frac{1}{\mu^{\prime}}\frac{\partial A^{\prime}}{\partial z}+M\cos \alpha_m\right]\\
  +&\frac{\partial}{\partial r} \left[ \frac{1}{\mu^{\prime}}\frac{\partial A^{\prime}}{\partial r}-M\sin \alpha_m\right]=-J
  \end{aligned}
\end{equation}
其中， $\mu^{\prime}=r\mu_0$， $A^{\prime}=rA$， $B^{\prime}=rB$， $\alpha_m$ 是 $\boldsymbol{M}$ 与 $r$轴正方向的夹角。 公式 \eqref{huanjian1}与二维平面情形非常相似，因此可以使用相同的方法进行处理 \cite{jinjianming}。
\subsubsection{边界条件}
%In Fig. \ref{bd}, the boundary condition is applied. The magnetic potential $A$ is zero in the infinite space, assuming that boundary dotted line is at infinity, therefore we have the first type boundary condition
在图 \ref{bd} 中，施加了相应的边界条件。在无限远的空间，磁势$A$为0，假定虚线边界为无限远，我们施加第一类边界条件
\begin{equation}\label{bd1}
  \Gamma_1:\,\,\,\,  A^{\prime}=0.
\end{equation}
%In addition, the $z$ axis is also applied with the second type boundary condition
除此之外，在$z$轴上，我们施加了第二类边界条件
\begin{equation}\label{bd2}
  \Gamma_2:\,\,\,\,  \frac{\partial A^{\prime}}{\partial \boldsymbol{n}}=0.
\end{equation}
\subsubsection{形函数}
%The Galerkin method \cite{polycarpou2005introduction} is employed to convert a differential equation to a discrete problem. %In Fig. \ref{ele},
伽辽金方法\cite{polycarpou2005introduction}用来将一个微分方程转化为一个离散问题来求解。
%Assuming $\triangle DEF$ is one of the triangles in the mesh network, in the heart of the triangle, the magnetic potential $A^{\prime}(r,z)$ is determined by the linear interpolation equation:
假定三角形$\triangle DEF$是有限元分网中的一个三角形，在三角形的内部，磁势$A^{\prime}(r,z)$由一下插值函数确定：
%\begin{figure}[t]
%  \centering{\includegraphics[width=4 cm]{ele}}
%  \caption{Triangular element in FEM}
%  \label{ele}
%\end{figure}
\begin{equation}\label{linear}
\begin{aligned}
  A^{\prime}(r,z)=\sum_i N_i(r,z)A^{\prime}_i \,\,\,\,(i=D,E,F),
  \end{aligned}
\end{equation}
其中，
\begin{equation}\label{Ni}
\begin{aligned}
  N_i(r,z)=\frac{1}{2\Delta} (a_i+b_ir+c_iz)\qquad(i=D,E,F)
  \end{aligned}
\end{equation}
是形函数，$\Delta$是三角形的面积。
\begin{equation*}\label{abc}
\left\{
\begin{aligned}
a_{ D}&=r_{ E}z_{ F}-z_{ E}r_{ F};\!\!\!&b_{ D}&=z_{ E}-z_{ F};\!\!\!&c_{ D}&=r_{ F}-r_{ E},\\
a_{ E}&=r_{ F}z_{ D}-z_{ F}r_{ D};\!\!\!&b_{ E}&=z_{ F}-z_{ D};\!\!\!&c_{ E}&=r_{ D}-r_{ F},\\
a_{ F}&=r_{ D}z_{ E}-z_{ D}r_{ E};\!\!\!&b_{ F}&=z_{ D}-z_{ E};\!\!\!&c_{ F}&=r_{ E}-r_{ D}.\\
\end{aligned}
\right.
\end{equation*}
%In the triangle element, the magnetic flux density in the direction of $r$ and $z$ is respectively written as
%\begin{equation}\label{Bxy}
%\left\{
%\begin{aligned}
%B_r=\frac{1}{2\Delta}\sum_i c_i A_i \\
%B_z=\frac{1}{2\Delta}\sum_i b_i A_i
%\end{aligned}
%\right.
%\quad(i=K,M,N).
%\end{equation}
\subsubsection{单元矩阵}
%Applying Galerkin method to \eqref{huanjian1}, we obtain the left side item where the unit coefficient matrix is
应用伽辽金方法到公式\eqref{huanjian1}，我们可以获得左侧的单元系数矩阵为
\begin{equation}\label{Se}
  S^eA^{\prime e}=
\left[
\begin{aligned}
S^e_{ DD} && S^e_{ DE} && S^e_{ DF}\\
S^e_{ ED} && S^e_{ EE} && S^e_{ EF}\\
S^e_{ FD} && S^e_{ FE} && S^e_{ FF}
\end{aligned}
\right]
\left[
\begin{aligned}
A^{\prime}_{ D} \\
A^{\prime}_{ E} \\
A^{\prime}_{ F}
\end{aligned}
\right]
\end{equation}
其中，
\begin{equation}\label{Sij}
\begin{aligned}
  S^e_{ij}=\frac{1}{4\Delta \mu^{\prime}}(c_ic_j+b_ib_j)\quad(i,j=D,E,F).\\
  \end{aligned}
\end{equation}
右侧电源项为
\begin{equation}\label{jm}
\begin{aligned}
J^e+M^e
=\frac{1}{3}J\Delta
\left[
\begin{aligned}
1 \\
1 \\
1
\end{aligned}
\right]+
\frac{M}{2}\left[
\begin{aligned}
c_{ D} \cos \alpha_m - b_{ D} \sin \alpha_m \\
c_{ E} \cos \alpha_m - b_{ E} \sin \alpha_m \\
c_{ F} \cos \alpha_m - b_{ F} \sin \alpha_m
\end{aligned}
\right].
\end{aligned}
\end{equation}
当我们把所有的单元装配一起，我们得到一个非线性的有限元方程：
\begin{equation}\label{femm}
  \boldsymbol{S}\boldsymbol{A^{\prime}}=\boldsymbol{J}+\boldsymbol{M}.
\end{equation}
\subsection{优化的TLM方法}
\subsubsection{等效电路模型}
%The idea of building an equivalent circuit is based on the fact that FEM's global matrix is obtained by assembling all elements' unit matrices. Once each unit matrix is modelled as a small circuit segment, we can group all these small pieces together, which means the FEM problem can be solved in a circuit way. In \eqref{jm}, the vector $J^e$ and $M^e$ can be easily modelled as current sources flowing into the nodes. For the unit coefficient matrix $S^e$ obtained in \eqref{Se}, it should be treated as a circuit nodal admittance matrix in TLM. It includes the admittances among the three nodes $D,E,F$. For example, for the node $D$, its self-admittance is
通过建立等效电路来进行有限元求解的想法来自于一个事实，那就是有限元的全局矩阵是通过将所有的单元矩阵进行装配得到的。如果每一个单元矩阵都被建模为一个小的电路单元，我们就可以将所有的电路单元组装起来，得到一个完整的电路网络。简言之，我们就能通过求解一个电路的方法来求解有限元问题。在公式\eqref{jm}当中，向量$J^e$和$M^e$可以很容易地建模为电流流入节点的电流源。对于公式\eqref{Se}中的单元系数矩阵$S^e$，在这里我们将它看做电路的节点导纳矩阵。它包含三个节点$D,E,F$之间的导纳。例如，对于节点$D$来说，它的自导为
\begin{equation}\label{ykk}
  Y^{e}_{ DD}=S^{e}_{ DD},
\end{equation}
节点$E$与节点$F$之间的互导为
\begin{equation}\label{ykm}
  \begin{matrix}
  Y^e_{ DE}=-S^e_{ DE}=-S^e_{ ED}, \\
  Y^e_{ DF}=-S^e_{ DF}=-S^e_{ FD}.
  \end{matrix}
\end{equation}
%To build a correct equivalent circuit model, we need to judge the sign of the admittances. First, it is necessary to check if there is other branch between node $D$ and the ground. Since
为了建立一个正确的等效电路模型，我们需要判断每一个导纳的正负号。首先，有必要检查一下存不存在节点$D$与地之间的分支。因为
\begin{equation}\label{ground}
  S^e_{ DD}+S^e_{ DE}+S^e_{ DF}=0,
\end{equation}
%the mutual admittance between node $D$ and the ground is zero, i.e., there is no branch connecting node $D$ and ground.
节点$D$与地之间的导纳为0，说明节点$D$与地之间没有分支。

%Then, we need to determine if the mutual admittances are positive or negative. Fortunately, it can be easily proven that
然后，我们需要确定导纳的正负。很容易证明
\begin{equation}
\left\{
\begin{aligned}
S^e_{ DE}=&\frac{1}{4\Delta \mu^{\prime}}(c_{ D}c_{ E}+b_{ D}b_{ E})=\frac{1}{4\Delta \mu^{\prime}}\overrightarrow{EF}\cdot\overrightarrow{DF}, \\
S^e_{ DF}=&\frac{1}{4\Delta \mu^{\prime}}(c_{ D}c_{ F}+b_{ D}b_{ F})=\frac{1}{4\Delta \mu^{\prime}}\overrightarrow{FE}\cdot\overrightarrow{ED}.
\end{aligned}
\right.
\end{equation}
%Evidently, the sign of $S^e_{DE}$ and $S^e_{DF}$ are separately related to $\angle DFE$ and $\angle DEF$, namely the angles that edge $DE$ and edge $DF$ faces respectively. As shown in Fig. \ref{triangle}, if the corresponding angle is less than $90^{\circ}$, the mutual admittance is positive, which should be modelled as a resistance. If the corresponding angle is greater than $90^{\circ}$, the mutual admittance is negative, which means that the equivalent component is an active component, then a voltage controlled current source is derived. If the corresponding angle is equal to $90^{\circ}$, there is no connection between the two nodes.
可见，$S^e_{DE}$与$S^e_{DF}$的正负与三角单元中$\angle DFE$ 和$\angle DEF$大小有关，即边$DE$与边$DF$所对的内角。如图\ref{triangle}所示，如果所对应的角度小于$90^{\circ}$，对应的互导为正，建模为纯电阻；如果所对应的角度大于$90^{\circ}$，对应的互导为负，该元件应为含源元件，由于它两端的电压是外部提供的，所以建模为电压控制电流源；如果所对应的角度等于$90^{\circ}$，对应的互导为零，表示两个节点之间没有电路连接。

%From the above analysis, we conclude that the pure resistance network only contains triangles in Fig. \ref{ruijiao} and Fig. \ref{zhijiao}. However, even if an efficient mesh program were used, the generated mesh network could inevitably contain obtuse angled triangle in Fig. \ref{dunjiao}. Of course, we could refine the mesh or redo the mesh, in which case a lot of time will be used into checking the mesh quality and modifying the mesh. In addition, for some small size elements, they can be difficult to refine. To sum up, the original resistances network should be generalized as a hybrid network of resistances $\boldsymbol{R}$, voltage controlled current sources $\boldsymbol{I_{c}}$ and current sources $\boldsymbol{J}$ and $\boldsymbol{M}$ for any shape triangles, and the circuit can be represented by the nodal admittance matrix
从上面的分析可以看出，一个纯电阻电路只包含图 \ref{ruijiao}和图 \ref{zhijiao}中的三角形。然而，即使是再优秀的分网程序，所生成的有限元网络也不可避免的出现钝角三角形，如图 \ref{dunjiao}所示。当然，我们可以通过优化分网或者重新进行分网来解决这个问题。但是这样就会将一些时间用于检查分网质量和修改分网上了。另外，对于一些非常小的单元，它们可能很难再被优化。所以，以前的研究学者所提出的纯电阻的等效电路应该被一般化为一个纯电阻$\boldsymbol{R}$和电压控制电流源$\boldsymbol{I_{c}}$和电流源$\boldsymbol{J}$、$\boldsymbol{M}$组成的混合网络，这个模型适合一般有限元三角分网单元。等效电路可以表达为
\begin{equation}\label{nam}
  \boldsymbol{Y}\boldsymbol{V}=\boldsymbol{J}+\boldsymbol{M}+\boldsymbol{I_{c}}
\end{equation}
%where $\boldsymbol{Y}$ is the matrix of admittances including unknown nonlinear components, $\boldsymbol{V}$ is node voltage matrix matching magnetic potential $A^{\prime}$ in \eqref{femm} and the right side is the current sources.
其中，$\boldsymbol{Y}$是电路的导纳矩阵，它包含非线性元件，$\boldsymbol{V}$是节点电压矩阵，跟公式 \eqref{femm} 中的磁势 $A^{\prime}$ 相对应。等侧右侧是电流源。

% In \cite{486548}, J. Lobry has mentioned that the reflection coefficient in pure resistance's TLM iteration is
%\begin{equation}\label{ref}
%  K_R=\frac{R-Z_0}{R+Z_0},
%\end{equation}
%where $Z_0$ is the transmission line's impedance. Normally, for passive components, $R>0$. However, for the current source, we have $R<0$. If \eqref{ref} is still used, the coefficient $K_R$ tends to be an extremely large number which can easily cause numerical instability when $Z_0$ is chosen as closely as the true value of current source's $R$ (in fact it's a common way to improve TLM's convergence speed). The FEM is an approximation solution of the real problem, in which the quality of mesh elements affects the solution error. In the mesh process, we should avoid generate elements with sharp angle, flat shape, or distorted shape. It's no doubt that
\begin{figure*}[t]
    \centering{
        \subfigure[锐角三角形]{
            \label{ruijiao} %% label for first subfigure
            \centering{\includegraphics[width=1.7cm]{ruijiao}}}
        \hspace{3mm}
        \subfigure[直角三角形]{
            \label{zhijiao} %% label for second subfigure
            \centering{\includegraphics[width=1.7cm]{zhijiao}}}
        \hspace{3mm}
        \subfigure[钝角三角形]{
            \label{dunjiao} %% label for second subfigure
            \centering{\includegraphics[width=2cm]{dunjiao}}}
            }
    \caption{有限元分网当中三种不同形状的的三角形
    }
    \label{triangle} %% label for entire figure
\end{figure*}
\begin{figure*}[t]
    \centering{
        \subfigure[非线性电阻的传输线模型]{
            \label{systemr} %% label for first subfigure
            \centering{\includegraphics[width=6cm]{systemr}}}
        %\hspace{2mm}

        \subfigure[入射阶段的等效电路]{
            \label{systeminduce} %% label for second subfigure
            \centering{\includegraphics[width=4.cm]{systeminduce}}}
        \subfigure[反射阶段的等效电路]{
            \label{systemreflect} %% label for second subfigure
            \centering{\includegraphics[width=3.5cm]{systemreflect}}}
            }
    \caption{传输线在电路网络中的应用
    }
    \label{sysmodel} %% label for entire figure
\end{figure*}
\subsubsection{松弛方法}
%The TLM iteration basically includes two stages: incidence and reflection. The working principle has been explained clearly in \cite{486548}. By adding a transmission line segment, the nonlinear components and linear network are separated and nonlinear solution of the circuit network can be implemented. As shown in Fig. \ref{systemr}, $R$ is the nonlinear resistance. By means of transmission line, the resistance is connected to a linear circuit. %As shown in Fig. \ref{tlmodel},
%Assuming that the voltage signal is propagating between $a$ and $c$, at a moment, when the incidence voltage $V_i$ travels from $c$ to $a$, if an observer stands at $a$, the equivalent circuit is shown in Fig. \ref{systeminduce}. At next moment, when the voltage signal arrives at $a$ and reflects to $c$, the reflection voltage $V_r$ can be obtained by
传输线迭代法一般包括入射和反射两个阶段。具体的工作原理在文献中已经阐释的比较清楚了\cite{486548}。在电路当中添加传输线可以实现非线性电路的求解。。如图\ref{systemr}所示，$R$为非线性电阻，通过传输线，与线性电路网络相连，假设电压信号在节点$a$与$c$之间传播，电压$V_i$从$c$入射到$a$，此时，如果站在$a$点观察，则此时传输线的电压入射到节点$a$，等效电路如图\ref{systeminduce}所示。在下一时刻，当电压信号到达节点$a$并反射到节点$c$，反射电压$V_r$为
\begin{equation}\label{vi}
  V_r=V_{ab}-V_i
\end{equation}
%where $V_{ab}$ is the voltage difference between $a$ and $b$.  In every iteration, the linear system remains unchanged. The LU decomposition method can be used to solve the system:
其中，$V_{ab}$为节点$a$和$b$之间的电压差。在每一个迭代中，线性电路保持不变。通常，使用LU分解法来进行系统的求解：
\begin{equation}\label{yvi}
  \boldsymbol{Y_{linear}}\boldsymbol{V}=\boldsymbol{J}+\boldsymbol{M}+\boldsymbol{I_{c}}+2\boldsymbol{V_i}\boldsymbol{Y_0}
\end{equation}
%where $\boldsymbol{Y_{linear}}$ is admittance matrix grouped by linear system and transmission line's admittances and $2\boldsymbol{V_i}\boldsymbol{Y_0}$ represents the virtual current sources in the equivalent circuits. Then an LU factorization $Y_{linear}=LU$ can be done at the first step and reused by the next steps. Thus, the voltage $V_{ab}$ and reflection voltage $V_r$ can be obtained.
其中，矩阵$\boldsymbol{Y_{linear}}$是由线性电路与传输线导纳组成的导纳矩阵。$2\boldsymbol{V_i}\boldsymbol{Y_0}$代表等效电路当中的虚拟电流源。在求解的第一步，对矩阵进行LU分解$Y_{linear}=LU$，这样，分解的结果可以被之后的迭代重复利用。这样，我们就能够得到电压$V_{ab}$和反射电压$V_r$。

%If an observer stands at $c$, the equivalent circuit is shown in Fig. \ref{systemreflect}. When the voltage signal arrives at $c$ and reflects to $a$, the incidence voltage of next time is given as
如果观察者站在点$c$，等效电路如图\ref{systemreflect}所示。当电压信号到达点$c$并反射回点$a$，下一时刻的入射电压为
\begin{equation}\label{vr}
  V_i=V_{cd}-V_r
\end{equation}
%where $V_{cd}$ is the voltage difference between $c$ and $d$.
其中，$V_{cd}$为节点$c$与节点$d$之间的电压。
%\begin{figure}[t]
%    \centering{
%        \subfigure[]{
%            \label{tl} %% label for first subfigure
%            \centering{\includegraphics[width=4cm]{tl}}}
%        %\hspace{2mm}
%
%        \subfigure[]{
%            \label{induce} %% label for second subfigure
%            \centering{\includegraphics[width=4cm]{induce}}}
%        \subfigure[]{
%            \label{reflect} %% label for second subfigure
%            \centering{\includegraphics[width=4cm]{reflect}}}
%            }
%    \caption{The basic transmission line model
%    \figfooter{a}{Lossless transmission line}
%    \figfooter{b}{Equivalent circuit during incidence phase}
%    \figfooter{c}{Equivalent circuit during reflection phase}
%    }
%    \label{tlmodel} %% label for entire figure
%\end{figure}
%In original TLM iteration, the nonlinear resistance $R$ is determined by solving a group of nonlinear equations with the help of Newton-Raphson method in each element \cite{486548}. However, there is a easier way to deal with it. The nonlinear part of $R$ is the permeability $\mu$. Since the voltage $V_a$ ($V_c$) and $V_b$ ($V_d$) will become equal when the iteration converges, we can take advantage of voltages obtained from \eqref{yvi} to directly calculate the magnetic flux density $B$ in the element using
在初始的TLM迭代方法中，非线性电阻$R$是通过牛顿迭代法求解一组非线性方程组来得到的\cite{486548}。然而，有一个更加容易地方法去处理这个问题。电阻$R$的非线性部分为磁导率$\mu$。既然电压$V_a$ ($V_c$) 和 $V_b$ ($V_d$) 在迭代收敛时会达到相等，我们可以利用这个电压直接计算出磁感应强度$B$。计算公式为
\begin{equation}\label{bbrbz}
  B=\sqrt{B_r^2+B_z^2}
\end{equation}
其中，
\begin{equation}\label{Bxy}
\left\{
\begin{aligned}
%B_r=\frac{1}{2\Delta}\sum_i c_i A_i \\
%B_z=\frac{1}{2\Delta}\sum_i b_i A_i \\
B_{r}&=-\frac{\partial A}{\partial z}\!\!&\approx&- \sum \limits_{i} \frac {c_iA_i}{2\Delta} \\
B_{z}&=\frac{A}{r}+\frac{\partial A}{\partial r}\!\!&\approx& \frac{\sum \limits_{i}A_i}{\sum \limits_{i}r_i}+\sum \limits_{i} \frac {b_iA_i}{2\Delta}\
\end{aligned}
\right.
(i=K,M,N),
\end{equation}
%and find the permeability $\mu$ according to the material's B-H curve. Then no more nonlinear iteration is needed and it can be thought as a kind of relaxation method. Voltages $V_{cd}$ and $V_i$ can be obtained by solving the small linear circuit.
求出$B$之后，通过查找材料的B-H曲线，来确定磁导率$\mu$。在这个过程中，没有任何的非线性求解部分，这可以被认为是一种松弛方法。这样电压$V_{cd}$和$V_i$就可以在小的线性电路当中进行求解。
%has a nonlinear relationship with $V_{cd}$ which can be written as $R=f(V_{cd})$ and the small nonlinear circuit is solved by. Generally, if we determine $R$ as a function with node voltage, that is, $R=f(V_{c},V_{d})$, we can use an approximate relationship to directly calculate $R=f(V_{a},V_{b})$,


%Generally, the TLM has the characteristics of unconditional stability because a DC circuit should have a steady-state over time. However, in our calculation, it may not be well suited for rapidly varying nonlinear function. If the varying slope is too great near the true value, a small change will produce a larger change in turn and the iterative process may tend to diverge or oscillate. Thus, in this situation, it is necessary to reduce or limit the changes in iteration using a damping factor $\alpha$. The variable $A$ in step $k$ will be determined as follows:
一般来说，TLM迭代能够无条件收敛。因为对于直流电路，随着时间的变化，都会达到一个稳定状态。然而，在我们的计算当中，发现对于变化非常陡峭的非线性方程来说，一个非常小的变化可能会反过来产生一个很大的反馈，从而出现振荡和不稳定的现象。在这种情况下，有必要通过一个衰减系数$\alpha$来减小或者限制在迭代中的变化大小。迭代步$k$中的变量$A$被重新更新为：
\begin{equation}\label{al}
  A_k = A_{k-1} + \alpha(A_k-A_{k-1}).
\end{equation}

%\begin{figure*}[t]
%    \centering{
%        \subfigure[]{
%            \label{sys1} %% label for first subfigure
%            \centering{\includegraphics[width=7.8cm]{sys1}}}
%            \hspace{10mm}
%        \subfigure[]{
%            \label{sys2} %% label for second subfigure
%            \centering{\includegraphics[width=8cm]{sys2}}}
%            }
%    \caption{The transmission line in circuit network
%    \figfooter{a}{Nonlinear resistance with lossless transmission line}
%    \figfooter{b}{Equivalent circuit during incidence process}
%    \figfooter{c}{Equivalent circuit during reflection process}}
%    \label{sysmodel1} %% label for entire figure
%\end{figure*}
\begin{figure*}[t]
    \centering{
        \subfigure[原始矩阵]{
            \label{lvl} %% label for first subfigure
            \centering{\includegraphics[width=4cm]{lvl}}}
        %\hspace{2mm}
        \subfigure[按级别排序的矩阵]{
            \label{lvl2} %% label for second subfigure
            \centering{\includegraphics[width=4cm]{lvl2}}}
        \subfigure[有限元中一个大小为 $4965\, \times \, 4965$ 的按照级别排序的下三角矩阵]{
            \label{levelLU} %% label for second subfigure
            \centering{\includegraphics[width=5cm]{levelLU}}}
            }
    \caption{级别调度方法
 }
    \label{level} %% label for entire figure
\end{figure*}
\subsection{加速方法}
%The TLM method evidently reduces the solution time in every single iteration benefitted from reusing the result of LU factorization. In addition, according to characteristics of this method, we adopts two acceleration methods.
由于能够重复利用LU分解的结果，TLM方法能够很明显地提高单步迭代的计算效率。除此之外，考虑到这种方法的相关特点，我们采用两种方法来对它进行加速。
\subsubsection{预处理器}
%The admittance $Y_0=1/Z_0$ of transmission line is an important parameter which affects the convergence speed. The closer $Y_0$ is to true value $Y_{true}$, the faster it converges. To obtain an ideal initial value of $Y_0$ while avoiding excessive execution time, we take advantage of a relatively coarse mesh which can be quickly solved by N-R. The acquired results can be mapped into a finer mesh network according to their geometry location as an initial value, thereby improving the speed of iteration.
传输线导纳$Y_0=1/Z_0$是一个影响传输线迭代法收敛速度的重要参数。$Y_0$选取越接近真值$Y_{true}$，收敛速度就越快。为了获得一个理想的$Y_0$的初始值，避免过多的计算时间，我们先用牛顿迭代法计算一个相对比较粗糙的有限元分网。通过几何映射，将所得到的计算结果映射到更加精细的分网上去，并作为迭代的初始值，从而提高迭代速度。
\subsubsection{并行加速}
%The solution time in nonlinear FEM is not only determined by the iteration steps, but also the time cost in each single iteration. Hence, it's necessary to study its parallel algorithm. The calculation in each step of TLM generally includes incidence phase and reflection phase. Obviously, in reflection phase, the solution is done in every single element which is of course suitable to be parallelized.
非线性有限元的计算时间不仅决定于迭代步数，而且还跟每个迭代步的时间有关。因此，很有必要研究其中的并行化算法。TLM的每一个迭代步一般包括入射阶段和反射阶段。很明显，在反射阶段，非线性求解可以放到每一个单元中进行，这非常适合并行化。

%Nevertheless, in incidence phase, because the LU decomposition has been done in the first step, the real work in this phase is solution of lower triangular matrix $L$ and upper triangular matrix $U$  using forward and backward substitutions. The parallelization of solving triangular system is always a difficulty. Researchers have done many studies on it \cite{higham1995stability,chen2016parallel,hogg2013fast,hogg2012fast}. Here, we adopted a method called level schedule to solve the triangular system in parallel.
然而，在入射阶段，因为LU分解已经在第一步完成了，所以，在这一阶段的真正工作是分别使用前向和后向消去进行下三角矩阵L和上三角矩阵U的求解。三角矩阵的并行化求解一直是一个难题。许多研究者都对这个问题进行了研究\cite{higham1995stability,chen2016parallel,hogg2013fast,hogg2012fast}。这里，我们采用一种叫做界别调度的方法来并行地求解这个三角系统。

%To solve the triangular system
为了求解三角系统
\begin{equation}\label{lxb}
  Lx=b,
\end{equation}
%we firstly define the level $l(i)$ of unknown variable $x(i)$. In $i$th $(i=1,2,\dots,n)$ row, for all the columns $j$th $(j=1,2,\dots,i)$ that satisfy the condition $L_{ij}\neq 0$ where $L_{ij}$ is the $(i,j)$th entry in matrix $L$, we determine $l(i)$ as
我们首先计算每一个未知变量$x(i)$的级别$l(i)$。在第$i$ $(i=1,2,\dots,n)$行，对于所有的列$j$ $(j=1,2,\dots,i)$里，满足条件$L_{ij}\neq 0$，我们定义$l(i)$为
\begin{equation}\label{xlevel}
  l(i)=1+\max_{j}l(j).
\end{equation}
%which is described in Algorithm \ref{alg1}.
%\begin{algorithm}[h]\label{alg1}
%\SetAlgoLined
%\KwResult{$level(i)\;(i=1,2,\dots,n)$}
% \For{$i=1 \;\KwTo\; n$}
%{
%    $level(i)=0$\;
%    \For {$j=1\;\KwTo\;i$}
%    {
%       \If{$L(i,j)\neq 0$}
%      {
%          $level(i)=\max(level(i),level(j))+1$\;
%      }
%    }
%}
% \caption{Calculate the $level(i)$ of each unknown $x(i)$}
%\end{algorithm}

%After obtaining the level of each unknown, we reorder the rows and columns according to their levels, thus the triangle problem becomes
获得每一个未知变量的级别之后，根据级别的大小，我们对行和列上下左右变化，新的三角问题变为
\begin{equation}\label{lxb2}
  L^{\prime}x^{\prime}=b^{\prime},
\end{equation}
%and the unknowns in the same level can be determined simultaneously.
在同一个级别的未知量可以同时被求解出来。

%Fig. \ref{lvl} and Fig. \ref{lvl2} illustrate a $12\, \times \, 12$ triangular matrix ordered by its levels before and after. The numbers in boxes are the level of each unknown in column. In Fig. \ref{lvl2}, unknowns in a group $\{x_1\}$, $\{x_2,x_5\}$, $\{x_3,x_6,x_9\}$, $\{x_4,x_7,x_{10}\}$, $\{x_8,x_{11}\}$, $\{x_{12}\}$ can be solved simultaneously.
如图 \ref{lvl}和图 \ref{lvl2}，分别为对未知量的级别进行排序前后的矩阵，在图\ref{lvl2}中，变量$\{x_1\},\{x_2,x_5\},\{x_3,x_6,x_9\},\{x_4,x_7,x_{10}\},\{x_8,x_{11}\},\{x_{12}\}$，可以依次被同时求解出来。图\ref{levelLU}为对一个大小为$4965\, \times \, 4965$的稀疏矩阵使用级别调度法排序之后的矩阵，黑色小点为矩阵当中的非零位置，两条虚线之间的行为同一个级别，经计算，该矩阵的级别为52。在每一个级别内，我们可以对其进行并行求解。红色正方形区域为对角元素，可以直接根据$b$值确定未知量的值，特别地，当L矩阵为单位上三角矩阵的话，正方形区域内无需求解。正方形下面为绿色和白色交替区域，计算完成正方形区域后，按照行的方向并行处理每行元素，并更新相应的b值，如此直至求解结束。U矩阵的计算算法与L矩阵类似，只需改为上三角阵即可。

%\begin{figure}[t]
%  \centering{\includegraphics[width=5cm]{levelLU}}
%  \caption{A $4965\, \times \, 4965$ lower triangular matrix in FEM ordered by its level}\label{levelLU}
%\end{figure}
%Fig. \ref{levelLU} describes a sparse lower triangular matrix in FEM of size $4965\, \times \, 4965$ reordered by the levels. The little dots represent the non-zero entries in matrix and each level is separated by the dashed line. By calculating, the number of level of this matrix is 52. The red square zone is the diagonal entries, by which unknowns can be solved by $x^{\prime}(i)=b^{\prime}(i)/L^{\prime}(i,i)$. Especially, when L is a unit lower triangular matrix, the entries in the diagonal positions is exactly the solutions. Below the red square zones, there are rectangle zones. Using the unknowns obtained from the upper red zone, we can update the $b^{\prime}(i)$ in the row direction which can be run in parallel. Defining $nLevel$  as the number of levels, $levelStart(i)$ the column number of level $i$ starts in $L^{\prime}$ and $levelEnd(i)$ the column number of level $i$ ends in $L^{\prime}$. The solution process of triangular system is described in Algorithm \ref{alg2}. The solution of the upper triangular matrix U is similar to that of L. %Fig. \ref{flowchart} shows the process of the complete modified TLM method.
%\begin{algorithm}
%\SetAlgoLined
%\KwResult{$level(i)\;(i=1,2,\dots,n)$}
% \For{$i=1 \;\KwTo\; n$}
%{
%    $level(i)=0$\;
%    \For {$j=1\;\KwTo\;i$}
%    {
%       \If{$L(i,j)\neq 0$}
%      {
%          $level(i)=\max(level(i),level(j))+1$\;
%      }
%    }
%}
% \caption{Calculate the level of each unknown}
% \label{alg1}
%\end{algorithm}
\begin{algorithm*}
\SetAlgoLined
\KwResult{$x^{\prime}$}
 \For {$i=1\;\KwTo\;nLevel$}
 {
    \tcc{parallel square zone in Fig. \ref{levelLU}.}
    \For {$j=levelStart(i)\;\KwTo\;levelEnd(i)$}
    {
        $x^{\prime}(j)=b^{\prime}(j)/L^{\prime}(j,j)$\;
    }
    \tcc{parallel rectangle zone in Fig. \ref{levelLU}.}
    \For {$k=levelEnd(i)+1\;\KwTo\;n$}
    {
        \For {$m=levelStart(i)\;\KwTo\;levelEnd(i)$}
        {
            $b^{\prime}(k) = b^{\prime}(k) - L^{\prime}(k,m)x^{\prime}(m)$\;
        }
    }
 }
 \caption{并行化的下三角求解算法 $L^{\prime}x^{\prime}=b^{\prime}$}
 \label{alg2}
\end{algorithm*}


%\begin{figure}
%  \centering{\includegraphics[width=7cm]{flowchart}}
%  \caption{Flowchart of the proposed modified TLM method}\label{flowchart}
%\end{figure}
\subsection{电磁力计算}
%Magnetic force is an important quantity when we want to analyze the actuator's behaviors. There are more than one way to compute the magnetic force, such as Lorentz method, Maxwell stress tensor, and virtual work method \cite{757938}. Here, a method named Coulomb virtual work (CVW) method \cite{1062812,1063232,894410} is exploited for the ease to use.
在分析制动器的特性过程当中，电磁力是一个重要的变量。对于电磁力的计算，有许多种方法，例如，洛伦兹力法、麦克斯韦应力法、虚功法等\cite{757938}。这些方法各有优缺点。本文采用了一种叫做Coulomb虚功法的电磁力计算方法来进行电磁力的计算\cite{1062812,1063232,894410}。
%\begin{figure}
%  \centering
%  \includegraphics[width=3cm]{cvw}
%  \caption{CVW}\label{cvw}
%\end{figure}

